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Motivated by a recent experiment directly measuring the current-phase relation (CPR) in graphene under 
the influence of a superconducting proximity effect, we here study the temperature dependence of the CPR in 
ballistic graphene SNS Josephson junctions within the the self-consistent tight-binding Bogoliubov-de Gennes 
(BdG) formalism. By comparing these results with the standard Dirac-BdG method, where rigid boundary 
conditions are assumed at the S|N interfaces, we show on a crucial importance of both proximity effect and 
depairing by current for the CPR. The proximity effect grows with temperature and reduces the skewness of the 
CPR towards the harmonic result. In short junctions (L < f ) current depairing is also important and gives rise 
to a critical phase cj> c < tt/2 over a wide range of temperatures and doping levels. 

PACS numbers: 74.45,+c, 74.50.+r 



I. INTRODUCTION 

Ever since the discovery of graphene 1 there has been inter- 
est in merging the unique two dimensional massless relativis- 
tic Dirac fermion spectrum of graphene with superconductiv- 
ity. A few years ago, graphene superconductor-normal metal- 
superconductor (SNS) Josephson junctions were successfully 
fabricated!^ Here superconductivity is induced in the S re- 
gions of the graphene through proximity to superconducting 
contacts deposited on top of these regions, and, through the 
Josephson effect, a finite supercurrent can exist even in the 
N region of the graphene. One of the key properties of a 
SNS junction is the current-phase relation (CPR), and, very 
recently, the first direct measurement of the CPR in graphene 
SNS junctions appearedP 

Theoretically, ballistic graphene SNS junctions were stud- 
ied even before the first experimental results. Several 
novel features were discovered, such as specular Andreev 
reflection, 7 and a finite supercurrent in undoped graphene at 
zero temperature, despite the point-like Fermi surface!^ These 
results, along with most theoretical results on the Joseph- 
son effect in graphene SNS junctionsj^^lhave employed the 
Dirac-Bogoliubov-de Gennes (DBdG) formalism, where the 
standard BdG formulation is applied to the Dirac spectrum. 7 
To solve the resulting equations, rigid boundary conditions 
are assumed for the superconducting order parameter A, such 
that A takes on a fixed, non-zero value in S, whereas it is 
zero in N. However, such an approach explicitly ignores any 
processes in which A is reduced in the S regions, such as 
proximity effect or depairing by current (see e.g. Ref. lfl"3l 
for a review). Moreover, with the Josephson current in the 
DBdG approach usually calculated as the current carried by 
sub-gap Andreev bound states, this method is limited to junc- 
tion lengths L < £, where £ is the superconducting coher- 
ence length. Other, related methods have been employed to 
relax some of the constraints in the original works, but they 
still apply rigid boundary conditions for A! 14 " 16 lln fact, in or- 
der to not apply any boundary conditions on A at the S|N 
interface, a self-consistent treatment is needed. One such 



method is the self-consistent tight-binding (TB) BdG formal- 
ism, where one only assumes that the superconducting con- 
tacts induces a pairing potential into the graphene and then 
solve self-consistently for A in the whole SNS structureFEED 
Not only does this approach give an explicit calculation of the 
full proximity effect, including current depairing, but it also 
results in a Josephson current appropriately calcul ated f rom 
this proximity effect. Previous works by the author^SED em- 
ploying this formalism at zero temperature showed on some 
corrections to the CPR compared to the DBdG results, as well 
as deviations for the critical current as function of junction 
length. 

Now, with experimental data on the CPR in graphene at 
hand, it is of large interest to theoretically map out the tem- 
perature dependence of the CPR in graphene. Very recently, 
some theoretical temperature dependent CPR results appeared 
using the DBdG approach, 16 but otherwise all theoretical work 
have been at zero temperature. The goals of this work are thus 
twofold: (1): Establish the correct CPR in ballistic graphene 
SNS Josephson junctions as function of temperature using the 
self-consistent TB BdG formalism. (2): Determine how well 
the DBdG approach captures the CPR, and, if applicable, de- 
termine the source(s) of discrepancy. 

We will here show that both the proximity effect depletion 
of superconductivity in the S regions and depairing by current 
are large in short junctions (L < £). In fact, these effects 
lead to a CPR where the critical current is reached for a phase 
4> c < it/2. This is true over a large temperature range and 
even more prominent for high doping levels in the graphene. 
In order to capture these effects a self-consistent solution of 
A is crucial since rigid boundary conditions explicitly ignores 
any such processes in the S regions. For longer junctions 
(L > £) the proximity effect depletion is unchanged but de- 
pairing by current is now not a big issue, and here <f> c > tt/2. 
In this junction length regime the DBdG approach using An- 
dreev bound states is not formally justified, and due to a large 
proximity effect, a self-consistent approach is still needed. 

The rest of the article is organized as follows. In the 
next section we briefly introduce the methods used, both the 
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self-consistent approach and the DBdG equation using rigid 
boundary conditions. Then in Section III we report our re- 
sults, including the temperature dependent CPR for different 
junctions as well as an analysis of the influence of proximity 
effect and current departing on the CPR. Finally, we summa- 
rize our key findings and comment on the applicability of our 
results to current and future experiments. 



II. METHOD 

Our starting point is the nearest neighbor tight-binding 
Hamiltonian on the graphene lattice together with an on-site s- 
wave superconducting order parameter Aj, which is induced 
by the proximity to external superconducting contacts: 

H = ~ * E ( a L^V + Hc ) + E Via(al a a ia + b\ a b ia ) 

i 

Here a\ a (b ia ) creates an electron with spin a on the A (B) 
sublattice in unit cell i. t ~ 2.5 eV is the nearest neighbor 
hopping amplitude, whereas is the chemical potential. We 
will assume that, due to the external contacts, the chemical 
potential is rather large and constant in the S regions whereas 
it can be tuned (to a constant value) with a back gate voltage 
in the N region. 



A. DBdG equation 

For an analytical treatment within the DBdG framework we 
set A N = 0, whereas A s = A(T) with A(T = 0) = A . 
Here A(T) is set to the standard BCS temperature depen- 
dence. This amounts to applying rigid boundary conditions 
on A, at the S|N interface. By treating the S and N regions 
separately, the above Hamiltonian can be Fourier transformed 
and the kinetic energy linearized around the Dirac points to 
produce the standard DBdG equation!^ 



H -fi A 
At -H + fJ, 



(2) 



where A and /i take on different, but constant, values in 
the S and N regions, as described above. Moreover, Ho = 
—ihvp(a x d x + <Jydy), where Hvp — V3ta/2 with a — 
2.46 A being the lattice constant of graphene. 

The strategy for calculating the Josephson current in the 
junction is to first obtain the energy spectrum for the Andreev- 
bound states in the N region. This is done by matching the 
wave functions given by Eq. ^ at the two S|N interfaces and 
then solve for the allowed energy eigenstates, e. The general 
expression for the Josephson current can then be written as! 2 -^ 



± n 



(3) 



where the summation over ± denotes all Andreev levels and 
the summation over n denotes the transverse momentum in- 
dex. Here, cf> is the phase drop across the junction,/(x) is 
the Fermi-Dirac distribution function, and the prefactor of 4 
is due to the spin- and valley-degeneracy. Let us first consider 
the case of /is = /ijv, i.e. no Fermi level mismatch (FLM) 
in the junction. In the wide junction limit where the junction 
width W satisfies W 3> L, we may replace the summation 
over discrete transverse momentum indices with an integral 
as follows: 



E 



2ir/W ~ 2kHv f 



d9 cos 9, 



(4) 



where the integration over the angle 9 takes into account all 
possible trajectories. The Andreev levels in a SNS junction 
with no FLM have the formP^ 



e± = ±A(T)y 1 - rsin 2 (0/2) 
where for graphene r has the specific form: 

cos 2 9 



1 — sin 9 cos 2 [/iatL cos 9/{hvp)] 



(5) 



(6) 



Inserting the above equations in to Eq. Q and introducing 
the normalization constant Io — (W/£)eAo/h, where £ = 
hvp/Ao is the superconducting coherence length, we arrive 
at the expression: 



fx n l 7 * 12 J/1 [A(T)] 2 rsin^cosf 



A9 



Io 27rAo7„ 7r/2 A e+ tanbT 1 ( ( 3£ + /2) ' 



(7) 



where ft is the inverse temperature. Finally turning our atten- 
tion to the FLM case, a similar procedure 8 leads to the expres- 
sion: 



I{4>) A(T)£ 



E 



(8) 



i a w ^ r, . 2,,/^ 

n y 1 - r n sin (0/2) 
where r n is given by Eq. (17) in Ref. (HJ. 

B. Self-consistent treatment 

In the self-consistent numerical treatment we do not prede- 
termine the value of A^ nor apply rigid boundary conditions 
at the S|N interface. Instead we assume that the influence of 
the superconducting contacts on the underlying graphene is 
only through an induced attractive pairing potential. For s- 
wave pairing the simplest potential is a constant, non-zero, 
attractive Hubbard-[/ term in the S regions. By applying the 
mean-field approximation to the attractive Hubbard model we 
reproduce Eq. ([T]) but with an added self-consistency condi- 
tion 



-Ui 



(audit) + (bi\b^) 



(9) 
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where Ui — U in S but zero in N. We can now solve self- 
consistently for Aj by first guessing a profile for A; through- 
out the SNS junction, solving Eq. ([TJi with this guess, recal- 
culating Aj using Eq. d9), and then reiterate this process until 
two subsequent A, profiles are within a predetermined error 
margin. In order to study the proximity effect between the S 
and N regions in the graphene we focus on the pairing am- 
plitude Fi = Ai/U{, While A, will only be non-zero in the 
S regions, Fi can, due to proximity effect leakage from S to 
N, be non-zero even in N. Alongside this leakage also comes 
a depletion of F in S such that F < A /f7 close to the in- 
terface. It is this latter process that, as we will show below, 
significantly changes the CPR. The other important property 
is the Josephson current I that flows through the junction if 
there exists a phase difference (f> in Aj across N. From the 
continuity equation for the charge current we can, using the 
self-consistent solution for Aj, calculate /. For additional de- 
tails on this self-consistent method see Ref. 1 19 1. 



C. Simulation details 

We will assume clean, smooth interfaces such that a Fourier 
transform along the direction parallel to the interface is appli- 
cable. In the self-consistent treatment, the type of interface 
can be varied but with an on-site pairing the direction of the 
interface will not matter. We will here use the zigzag inter- 
face, such that one unit cell is \fZaj2 long. Moreover, we will 
only consider the wide junction regime, W 3> L ~ £. In the 
self-consistent treatment the junctions are naturally infinitely 
wide, due to Fourier transforming along the interface, whereas 
in the DBdG approach W = 30£ is used. For a direct com- 
parison between the two methods the current will be given in 
units of W/C 

In this article we will focus on a few representative val- 
ues for the various physical input parameters. In the super- 
conducting S regions we use U = 1.36t and /is = 0.6t. 
These values give Ao = 0.042i, and a superconducting co- 
herence length £ = hvpj Ao ~ 24 unit cells. This satisfies 
£ <C Xp — fivp/n in both S and N, a requirement for the 
DBdG solution, and will also allow us to self-consistently in- 
vestigate both the L < £ and L > £ regimes, where L is the 
length of the N region. Moreover, with these values we get 
S regions as small as 50 unit cells displaying clear bulk be- 
havior, which is advantageous due to the the computational 
demands of the self-consistent method. Note though that for 
a direct comparison with an experimental setup A is very 
large. However, we believe that our results are applicable even 
for smaller Aol^and we, therefore, explicitly only report our 
results in dimensionless units. We will vary the chemical po- 
tential in the N region from the Dirac point (/Ltjy = 0) to mod- 
erately doped (hn — 0.47/15), to no FLM at the interface 
(/j at = /is). In addition, we will investigate both short junc- 
tions with L = 0.42£, where our DBdG approach is techni- 
cally justified, and long junctions with L = 2.1£. 

Finally, let us comment on the application of a phase gra- 
dient in Aj. We apply a phase difference A<fi between the 
superconducting order parameters in the outermost regions of 



the two S regions. In the self-consistent treatment the phase 
is, however, allowed to relax also in the S regions close to the 
interfaces. Thus the phase drop <\> across the junction itself, 
i.e. across N, will necessarily satisfy <p < A<fr < tt in the 
self-consistent treatment, with an inequality in the first step as 
soon as the supercurrent is non-zero. As a consequence, we 
will not be able to trace out the self-consistent CPR for <f> close 
to n injunctions with high currents. While this appears as a 
numerical artifact in this context, it is, in fact, closely related 
to the physical 2ir phase-slip process in Josephson junctions 
(see, e.g., Ref. ll22lH . 



HI. RESULTS 

We will start by reporting the CPR for several different bal- 
listic graphene SNS junctions. Then we will explicitly ana- 
lyze the critical current and phase, followed by a qualitative 
analysis of the effect of proximity effect and current depair- 
ing on the CPR. The latter two quantities are captured by our 
self-consistent approach but is not included within the DBdG 
framework. 



A. CPR 

Figure [T] shows the CPR for four representative cases with 
/j at = and no FLM, and for short and long junctions, re- 
spectively. For each case we plot the CPR at four different 
temperatures, T/T c = 0, 0.14, 0.43, and 0.87. For short junc- 
tions we also report the DBdG results as dashed curves for a 
straightforward comparison. We will comment more specif- 
ically on the skewness, or anharmonicity, of the CPR when 
extracting cf> c in Fig. [2] but we see directly that the DBdG re- 
sults only approximately reproduces the self-consistent results 
in the limit of low temperature and low doping levels in N. 
For increasing doping levels and/or increasing temperatures, 
the self-consistent results have a skewness towards C < ir/2, 
whereas 4> c > ir/2 is always the case for the DBdG results. 
By parameterizing the skewness as S — 2cj) c /ir — 1, we will 
refer to the former case as negative skewness. We also directly 
see that the total current in the junction is naturally increased 
with increasing doping level, decreasing junction length, or 
decreasing temperature. Note though that the self-consistent 
results show a much larger suppression of the current with in- 
creasing temperature than the DBdG results. 

Before we continue, the case of no FLM at the interface 
deserves special attention. Here there is no interface barrier 
and we thus have an SNS junction with fully transparent in- 
terfaces. Such junctions, independent on junction length, have 
been shown to have a 1/T diverging superconducting decay 
length £at in the normal region. 18 For T = this is directly 
manifested in a linear drop of A<f) throughout the whole junc- 
tion (i.e. even in the S regions) and, consequently, also a linear 
CPR, resulting in <f> c — > tt. 19 The diverging £y at T = also 
leads to currents not limited by tunneling through Andreev 
bound states but by the size of the (induced) superconducting 
gap in the N region. Due to these high currents we cannot 
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FIG. 2: (Color online) Critical phase (f> c /n (a, b) and critical cur- 
rent 7c /To (c, d) as a function of the reduced temperature T/T c for 
/us = (black, x), 0.47 (green, o), and no FLM (red, □) for 
short junctions (a, c) and long junctions (b, d). Self-consistent results 
(solid, lines are only guides to the eye), and DBdG results (dashed). 
Note the log scale in (c, d). 



FIG. 1: (Color online) Josephson current I in units of To = 
(W/(,)eAo/h, when vf = 1, as a function of the phase drop (f> 
across the junction for /ijv = (a, b) and no FLM (c, d) for short 
junctions (a, c) and long junctions (b, d) at T/T c = (black, x), 
0.14 (blue, o), 0.43 (red, □), and 0.87 (green o). In (d) T/T c = 
0.025 (magenta, A) is also plotted whereas it falls on top of the 
T = results in the other plots. Self-consistent TB BdG results 
(solid, lines are only guides to the eye), and DBdG results (dashed). 



fully trace out the CPR in these junctions at low tempera- 
tures. However, in the accessible range, the self-consistent 
CPR is indeed linear, which supports the claim that <fi c — > ir as 
T — >• 0. Note that this numerical problem can partially be cir- 
cumvented by studying longer junctions, as seen in Fig.[TJd), 
since these junctions carry smaller currents due to a smaller 
superconducting gap in the N region. Since the current cal- 
culation in the DBdG framework explicitly assumes the pres- 
ence of a finite number of Andreev bound states in the normal 
(metal) N region, this method cannot capture the effects of the 
diverging £/v". This results in an underestimation of both the 
skewness and the absolute value of the critical current injunc- 
tions with no FLM at very low temperatures when using the 
DBdG method. 

Let us now focus on the critical phase <j> c , extracted for 
three different doping levels in Figs. [2] (a, b). In short junc- 
tions the DBdG results (dashed lines) show a notable (posi- 
tive) skewness away from the traditional harmonic Josephson 
relation I = I c s'm((j>). At the Dirac point at zero temperature, 
4> c = 0.637T as already established analytically, 8 but with in- 
creasing temperature the CPR approaches the harmonic form, 
in agreement with other recent DBdG results. 16 The skewness 
in the DBdG approach is somewhat increased with increasing 
doping but the overall temperature dependence is not highly 
sensitive to the doping level in N, nor is the no FLM case dis- 
tinctly different from the other results. What is most striking 
with Fig. |2|a) is the fact that the DBdG results do not repro- 



duce the self-consistent results reliably. At low doping levels, 
4> c is essentially uniformly shifted to lower values in the self- 
consistent approach, such that (f> c < tt/2 above some temper- 
ature. The temperature at which <j) c = tt/2 decreases with 
increasing doping levels. This leads to a substantial negative 
skewness with increasing doping, and thus the self-consistent 
results deviate even more from the DBdG results as the doping 
increases in the N region. At very high temperatures there is 
finally a tendency to approach a harmonic CPR, which is the 
development for any ballistic SNS junction with rigid bound- 
ary conditions as T — > T C P Within our numerical accuracy 
we cannot, however, model temperatures close enough to T c 
to formally see if 4> c —> tt/2 as T — >• T c . 

For long junctions the skewness is instead in general pos- 
itive or only minutely negative for the highest doping levels. 
Also, <fi c — > tt/2 even for only moderately high temperatures. 
We can for L > £ not formally justify the use of our DBdG 
framework, but a comparison still shows on a large discrep- 
ancy between the self-consistent results and the DBdG results, 
although it is smaller than for shorter junctions. Notably, the 
self-consistent results reach the harmonic (f> c = tt/2 result 
much faster than the DBdG results. 

In terms of the critical current, even with the log-scale in 
Figs.|2jc, d), we see that a self-consistent treatment results in 
a significantly lower current for all junctions except at very 
low doping levels and temperatures. In fact, the discrepancy 
between the two methods grow strongly with increasing tem- 
perature. 

In aggregate, the above results directly leads to two con- 
clusions: (1) a self-consistent approach is crucial in the short 
junction regime and still important for longer junctions. (2) 
a no FLM junction at zero temperature has a diverging su- 
perconducting decay length £at and is therefore fully super- 
conducting. This leads to a linear CPR with <p c — > tt as 
T —> 0. The large difference between the self-consistent and 
the DBdG results must stem from processes that are not cap- 
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tured when applying rigid boundary conditions to Aj at the 
S|N interface. We have already discussed the special case of 
no FLM junctions at zero temperature where A jy becomes 
finite. Also, very generally, processes influencing the super- 
conducting state in the S region are present in all junctions and 
at all temperatures. These processes include proximity effect 
depletion of the pairing amplitude F in the S region close to 
the interface and the additional loss of pairs in the whole junc- 
tion because of a finite supercurrent, or, equivalently, a finite 
phase gradient. In fact, shifts of <fi c to the region <fi c > ir/2 are 
mainly governed by processes in the N region, as evident in 
the DBdG approach, whereas shifts towards <fi c < ir/2 occur 
because of processes in the S regions. 13 That is, in the absence 
of the loss of Cooper pairs in the S regions of the interface, 
the self-consistent results would reproduce the DBdG results, 
where the intrinsic properties of the graphene N region lead to 
C > tt/2. 



B. Proximity effect 

Figure [3] contains the evidence for proximity effect deple- 
tion of pair amplitude at the interfaces. In (a) we see that for 
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FIG. 3: (Color online) (a): Absolute value of the pairing amplitude 
F in units of Fo = Ao/U at the critical current in a short junction 
for /ijv = (black) and no FLM (red) for temperatures T/T c = 
0, 0.43, and 0.87 (decreasing amplitude). Vertical lines mark the S|N 
interfaces. (The curve for no FLM at T = is at the highest achievable 
4>, which did not reach I c .) (b): Ratio of \F\ at the S|N interface to 
| F\ in the bulk as a function of the reduced temperature T/T c for jitjv 
= in a short junction (black, x) and long junction (blue, o) and for 
no FLM in a short junction (red, □) and long junction (green, 0). 

strong doping (red) the proximity effect is clearly larger, with 
more depletion of F in S and the accompanied accumulation 
of pairs in N. This is expected since with increasing doping 
in N, the effective barrier at the interface due to the FLM is 
decreased. The increase in the Josephson current with dop- 
ing of N is a direct consequence of this enhancement of F 
in N. In terms of temperature dependence, the pair amplitude 
reduction in the bulk S regions with increasing temperature 
follow the standard BCS temperature dependence and is thus 
also included in the DBdG treatment. On top of this, there 
is, however, an additional loss of pair amplitude both on the 



S and N sides of the S|N interfaces with increasing tempera- 
ture. This latter process is only captured in the self-consistent 
method. Also seen in (a) are oscillations in Fi at both ends 
of the S regions. These are due to the outer edges of the S 
regions, but do not influence the results of the junction itself. 
The small oscillations for [if? = at the interfaces at low 
temperatures are directly correlated with Friedel-like charge 
oscillations, which are only present when N is close to the 
Dirac pointful 

In (b) we analyze the proximity effect at the interface more 
closely. For strong doping in N (red, green), the proximity ef- 
fect is notably stronger than for an undoped N region (black, 
blue). However, note that the proximity effect is essentially 
independent on the junction length L. This directly tells us 
that the stronger discrepancy between the DBdG and the self- 
consistent results found for short junctions is not due to the 
proximity effect. Moreover, there is also a significantly in- 
creased interface proximity effect with increasing tempera- 
ture. This explains the stronger decrease in both the critical 
current and phase with increasing temperatures than found in 
the DBdG results. 



C. Current depairing 

Finally in Fig. [4] we study the absolute value of the pair 
amplitude for both = and C in short junctions. A 
difference in \F\ for these two phase gradients will directly 
signal the influence of a finite phase gradient or, by extension, 
a supercurrent on F. In the case of a loss of \F\ going from 
= and C this is known as depairing by current, a process 
present in any junction carrying a finite supercurrent. 
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FIG. 4: (Color online) Absolute values of the pairing amplitude F 
in units of Fo = Ao /U in short junctions for = (black) and 
(j> — <j) c (red) at two different temperatures for /in = (a) and no 
FLM (b). Vertical lines mark the S|N interfaces. 

In (a) we see that for low doping levels in N, there is a 
small amount of current depairing at low temperatures, which 
decreases with increasing temperatures. For heavily doped 
N regions (b), the depairing is significantly larger, though it 
also decreases with increasing temperature. The same set of 
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curves for long junctions show no change with applied phase 
gradient and thus depairing by current is only an issue in short 
junctions. This directly leads to the conclusion that it is cur- 
rent depairing which is responsible for the smaller and even 
negative skewness found in shorter junctions. Its doping de- 
pendence also explains why the skewness is more negative for 
larger doping levels. Moreover, the strong increase in proxim- 
ity effect with increasing temperature is partially counteracted 
in short junctions by the decreased current depairing. Thus we 
see less of a temperature dependence for the skewness when 
the junction length is decreased. In passing, we also note that 
the strength of the depairing process continues to increase as 
<f> is increased from <fi c towards tt, despite a decrease in the 
supercurrent. Thus the depairing is technically a function of 
the finite phase difference over the junction and not per see of 
the current, although they are intimately related. 



IV. CONCLUSIONS 

We have above self-consistently established the temper- 
ature dependent CPR in ballistic graphene SNS Josephson 
junctions with varying doping levels and junction lengths. 
We have also compared our results with the currently preva- 
lent method of studying these junctions, the so-called DBdG 
method, which explicitly relies on the use of rigid bound- 
ary conditions for the superconducting order parameter A at 
the S|N interfaces. The DBdG method gives a critical phase 
4> c = 0.6 — 0.7 at zero temperature, with larger values for 
higher doping levels in N, which slowly then reduces to the 
harmonic value of <fi c = it/2 as T — > T c . Our DBdG results 
are formally only valid in the short junction regime (L < £), 
but other recent work within the same DBdG framework has 
shown on a similar behavior for long junctions!^ The self- 
consistent results, however, are in most junctions qualitatively 
different from the DBdG results. First of all, any process influ- 
encing the superconducting state in the S regions of the junc- 
tion can only be captured when the rigid boundary conditions 
are removed and A is calculated self-consistently. Such pro- 
cesses are known 13 to decrease the skewness of the CPR and 
can even produce negative skewness, i.e. give <fi c < ir/2. We 
have studied both the influence of proximity effect depletion 
of A in S close to the interfaces and the additional loss of 
Cooper pairs when a finite phase gradient is applied across 
the junction. The latter process, is known as depairing by cur- 
rent and is present whenever there is a finite supercurrent in 
the junction. The proximity effect is important and modifies 
the CPR in all junctions we have studied. It increases with 
doping and, even more so, with temperature. It is, however, 
independent of the junction length. Current depairing, on the 
other hand, is only important in short junctions. It increases 
with doping, but decreases with temperature, both natural con- 
sequences of the amount of current in the junction. In fact, 
current depairing helps explain the growing discrepancy in I c 
with decreasing L between the analytical and self-consistent 
results found in Ref. |20|. Based on these results we can qual- 
itatively explain the modified CPR in the self-consistent re- 
sults: For long junctions only the proximity effect modifies the 



DBdG results. The proximity effect reduces the skewness sig- 
nificantly and also heavily suppresses the critical current, es- 
pecially at higher temperatures. However, for T approaching 
T c , cf) c is pinned to the harmonic result, thus implying that in 
the very high temperature limit rigid boundary conditions are 
still applicable. In shorter junctions, the same proximity effect 
is at play, but here current depairing causes further decrease 
in 4> c . This is especially true for higher doping levels where 
we find <f> c as small as 7r/4. The temperature dependence for 
the current depairing is opposite that of the proximity effect, 
and we thus see a somewhat smaller temperature variation in 
(f> c for shorter junctions. Also, for short junctions a noticeable 
anharmonicity is present closer to T c , which by itself strongly 
signals that these junctions are not well modelled using rigid 
boundary conditions. 13 In addition to these processes, junc- 
tions with no FLM at the interfaces have cj) c —> tt as T — )• 0, 
due to a 1/T diverging superconducting decay length <^v! 18 l 19 l 
Since the DBdG framework assumes Ajv = 0, also this effect 
is only captured in our self-consistent treatment. 

Let us finally also briefly comment on the applicability of 
our results to an experimental setup. In the recent exper- 
iment measuring the CPR in graphene SNS junctions'^ the 
skewness was found to always be positive and increasing lin- 
early with the critical current (i.e. with decreasing temperature 
and/or increasing doping). It was, however, concluded that 
the junctions were very likely in the quasi-diffusive regime 
and, thus, an explicit comparison with our work is not pos- 
sible. Nonetheless, in Ref. [16] a qualitative comparison was 
made to ballistic DBdG results in the appropriate long junc- 
tion regime (L ~ 3.50, which also showed linear behavior 
for small I c , although the slope did not match the experimen- 
tal results. Our results show that, even in this long junction 
limit, self-consistency is necessary in ballistic junctions in or- 
der to correctly capture the CPR, due to the proximity effect. 
The temperatures used in the experiment was, however, large 
enough that within our current model setup we cannot study 
these high temperatures, since A(T) becomes too small. Still, 
we can see that at slightly lower temperatures the skewness is 
not hugely temperature or doping dependent. We thus specu- 
late that when relaxing the ballistic requirement, proximity ef- 
fect corrections become less important. Ballistic transport has, 
however, already been demonstrated in graphene SNS junc- 
tions with the measurement of multiple Andreev reflections. 3 
Thus, a measurement of the CPR in ballistic graphene junc- 
tions does not seem to be a too distant an achievement. Our 
main, and easily verifiable, prediction in ballistic junctions - 
that the skewness becomes negative due to a combination of 
current depairing and proximity effect in short junctions - will 
therefore be important as the experimental junction lengths 
shrink in the future. 

In summary, we have found that proximity effect, depairing 
by current, and diverging superconducting decay length all 
qualitatively modifies the CPR in graphene Josephson junc- 
tions. Thus, a correct description of the CPR necessarily re- 
quires a self-consistent treatment of the superconducting or- 
der parameter. This is different from metallic SNS junctions 
where it has been shown that rigid boundary conditio ns ar e 
applicable for junctions with large FLM at the interface 1121221 
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